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Our preceding paper, Ref. 1, introduced a method to compute Casimir forces in arbitrary ge- 
ometries and for arbitrary materials that was based on a finite-difference time-domain (FDTD) 
scheme. In this manuscript, we focus on the efficient implementation of our method for geometries 
of practical interest and extend our previous proof-of-concept algorithm in one dimension to prob- 
lems in two and three dimensions, introducing a number of new optimizations. We consider Casimir 
piston-like problems with nonmonotonic and monotonic force dependence on sidewall separation, 
both for previously solved geometries to validate our method and also for new geometries involving 
magnetic sidewalls and/or cylindrical pistons. We include realistic dielectric materials to calculate 
the force between suspended silicon waveguides or on a suspended membrane with periodic grooves, 
also demonstrating the application of PML absorbing boundaries and/or periodic boundaries. In 
addition we apply this method to a realizable three-dimensional system in which a silica sphere is 
stably suspended in a fluid above an indented metallic substrate. More generally, the method allows 
off-the-shelf FDTD software, already supporting a wide variety of materials (including dielectric, 
magnetic, and even anisotropic materials) and boundary conditions, to be exploited for the Casimir 
problem. 



I. INTRODUCTION 

The Casimir force, arising due to quantum fluctua- 
tions of the electromagnetic field [2], has been widely 
studied over the past few decades [3j HJ El E] and ver- 
ified by many experiments. Until recently, most works 
on the subject had been restricted to simple geometries, 
such as parallel plates or similar approximations thereof. 
However, new theoretical methods capable of comput- 
ing the force in arbitrary geometries have already be- 
gun to explore the strong geometry dependence of the 
force and have demonstrated a number of interesting ef- 
fects [Il[8j|9l[10l[IIl[l2l[13l[T4l^ [16]. A substan- 
tial motivation for the study of this effect is due to re- 
cent progress in the field of nano-technology, especially 
in the fabrication of micro-electro-mechanical systems 
(MEMS), where Casimir forces have been observed [17] 
and may play a significant role in "stiction" and other 
phenomena involving small surface separations. Cur- 
rently, most work on Casimir forces is carried out by 
specialists in the field. In order to help open this field to 
other scientists and engineers, such as the MEMS com- 
munity, we believe it fruitful to frame the calculation of 
the force in a fashion that may be more accessible to 
broader audiences. 

In Ref. 1, with that goal in mind, we introduced 
a theoretical framework for computing Casimir forces 
via the standard finite-difference time-domain (FDTD) 
method of classical computational electromagnetism [18] 
(for which software is already widely available). The pur- 
pose of this manuscript is to describe how these compu- 
tations may be implemented in higher dimensions and 
to demonstrate the flexibility and strengths of this ap- 
proach. In particular, we demonstrate calculations of 
Casimir forces in two- and three-dimensional geometries, 
including three-dimensional geometries without any ro- 
tational or translational symmetry. Furthermore, we de- 



scribe a harmonic expansion technique that substantially 
increases the speed of the computation for many systems, 
allowing Casimir forces to be efficiently computed even 
on single computers, although parallel FDTD software is 
also common and greatly expands the range of accessible 
problems. 



Our manuscript is organized as follows: First, in 
Sec. [Il| we briefly describe the algorithm presented 
in Ref. [T| to compute Casimir forces in the time do- 
main. This is followed by an important modification 
involving a harmonic expansion technique that greatly 
reduces the computational cost of the method. Second, 



Sec. Ill presents a number of calculations in two- and 



three-dimensional geometries. In particular, Sec. Ill A 
presents calculations of the force in the piston-like struc- 
ture of Ref. O and these are checked against previous 
results. These calculations demonstrate both the va- 
lidity of our approach and the desirable properties of 
the harmonic expansion. In subsequent sections, we 
demonstrate computations exploiting various symmetries 
in three dimensions: translation-invariance, cylindrical 
symmetry, and periodic boundaries. These symmetries 
transform the calculation into the solution of a set of two- 



dimensional problems. Finally in Sec. HIE we demon- 
strate a fully three-dimensional computation involving 
the stable levitation of a sphere in a high-dielectric 
fluid above an indented metal surface. We exploit a 
freely available FDTD code [19 , which handles sym- 
metries and cylindrical coordinates and also is script- 
able/programmable in order to automatically run the se- 
quence of FDTD simulations required to determine the 
Casimir force [2D] . Finally, in the Appendix, we present 
details of the derivations of the harmonic expansion and 
an optimization of the computation of g(t). 
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FIG. 1: Schematic showing the two-dimensional piston- like 
configuration of [22] : Two perfectly conducting circular cylin- 
ders of radius R, separated by a distance a, are sandwiched 
between two perfectly conducting plates (the materials are 
either perfect metallic or perfect magnetic conductors). The 
separation between the blocks and the cylinder surface is de- 
noted as h. 



II. HARMONIC EXPANSION 

In this section we briefly summarize the method 
of Ref. [TJ and introduce an additional step which greatly 
reduces the computational cost of running simulations in 
higher dimensions. 

In Ref. 1, we described a method to calculate Casimir 
forces in the time domain. Our approach involves a mod- 
ification of the well-known stress-tensor method [21], in 
which the force on an object can be found by integrat- 
ing the Minkowski stress tensor around a surface S sur- 
rounding the object ( Fig. [I]), and over all frequencies. 
Our recent approach pQ abandons the frequency domain 
altogether in favor of a purely time-domain scheme in 
which the force on an object is computed via a series 
of independent FDTD calculations in which sources are 
placed at each point on S. The electromagnetic response 
to these sources is then integrated in time against a pre- 
determined function g(—t). 

The main purpose of this approach is to compute the 
effect of the entire frequency spectrum in a single simula- 
tion for each source, rather than a separate set of calcu- 
lations for each frequency as in most previous work [T2] . 

We exactly transform the problem into a mathemati- 
cally equivalent system in which an arbitrary dissipation 
is introduced. This dissipation will cause the electromag- 
netic response to converge rapidly, greatly reducing the 
simulation time. In particular, a frequency-independent, 
spatially uniform conductivity a is chosen so that the 
force will converge very rapidly as a function of simula- 
tion time. For all values of <r, the force will converge to 
the same value, but the optimal a results in the short- 
est simulation time and will depend on the system under 
consideration. Unless otherwise stated, for the simula- 
tions in this paper, we use a = 1 (in units of 27rc/a, a 
being a typical length scale in the problem). 

In particular, the Casimir force is given by: 



the electromagnetic fields on the surface S defined in our 
previous work pQ. 

Written in terms of the electric field response in direc- 
tion i at (£, x) to a source current J(£, x) = 5(t)5(x — x') 
in direction j, Eij{t\ x, x'), the quantity Tf(t) is defined 
as: 

Tf(t) EE jfd^x)^^ 

ee J ^(x)r^.(t;x,x) 

where dSj(x) = cLS (x) rij (x) , dS(x) the differential area 
element and n(x) is the unit normal vector to S at x. A 
similar definition holds for Tf(t) involving the magnetic 
field Green's function Hij. 

As described in Ref. [TJ computation of the Casimir 
force entails finding both T jE; (t;x, x) and the r^(t;x, x) 
field response with a separate time-domain simulation for 
every point x G S. 

While each individual simulation can be performed 
very efficiently on modern computers, the surface S will, 
in general, consist of hundreds or thousands of pixels or 
voxels. This requires a large number of time-domain sim- 
ulations, this number being highly dependent upon the 
resolution and shape of 5, making the computation po- 
tentially very costly in practice. 

We can dramatically reduce the number of required 
simulations by reformulating the force in terms of a har- 
monic expansion in T jE; (t;x, x), involving the distributed 
field responses to distributed currents. This is done as 
follows [an analogous derivation holds for r H (t;x, x)]: 

As S is assumed to be a compact surface, we can 
rewrite r^ (t;x, x) as an integral over S: 



rg(t;x,x)= f d5(x')r^(t;x,x')^(x-x') 

J S 



(2) 



where in this integral dS is a scalar unit of area, and 5s 
denotes a ^-function with respect to integrals over the 
surface S. Given a set of orthonormal basis functions 
{/ n (x)} defined on and complete over 5, we can make 
the following expansion of the S function, valid for all 
points x, x' G S: 



«5 s (x-x') = ^/„(x)/„(x') 



(3) 



The / n (x) can be an arbitrary set of functions, assuming 
that they are complete and orthonormal on S [31]. 

Inserting this expansion of the ^-function into Equa- 
tion ([2| and rearranging terms yields: 



Fi = Im- dt g(-t) (rf (t) +Tf (t)) , (1) r ^ x > x ) = (^^(x')r^(i;x,x')/„(x') 

7T Jq n 



where g(t) is a geometry-independent function discussed 
further in the Appendix, and the F(t) are functions of 



(4) 

The term in parentheses can be understood in a physical 
context: it is the electric-field response at position x and 
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FIG. 2: Differing harmonic expansions of the source currents 
(red) on the surface S. The left part shows an expansion us- 
ing point sources, where each dot represents a different sim- 
ulation. The right part corresponds to using / n (x) ~ cos(x) 
for each side of S. Either basis forms a complete basis for all 
functions in S. 



time t to a current source on the surface S of the form 
J(x, t) = S(t)f n {x). We denote this quantity by Vf-. n : 



r&„(t,x)= / ds(x' 



ir£(t;x,x')/„(x') 



(5) 



where the n subscript indicates that this is a field 
in response to a current source determined by / n (x). 
r^ ;n (£,x) is exactly what can be measured in an FDTD 
simulation using a current J(x, t) = 5(t)f n (x) for each n. 
This equivalence is illustrated in Fig. [2] 

The procedure is now only slightly modified from the 
one outlined in Ref. [T| after defining a geometry and 
a surface S of integration, one will additionally need to 
specify a set of harmonic basis functions {f n ( x )} on S. 
For each harmonic moment n, one inserts a current func- 
tion J(x,t) = S(t)f n {x) on S and measures the field re- 
sponse r^- ;n (x, t). Summing over all harmonic moments 
will yield the total force. 

In the following section, we take as our harmonic 
source basis the Fourier cosine series for each side of S 
considered separately, which provides a convenient and 
efficient basis for computation. We then illustrate its ap- 
plication to systems of perfect conductors and dielectrics 
in two and three dimensions. Three-dimensional sys- 
tems with cylindrical symmetry are treated separately, 
as the harmonic expansion (as derived in the Appendix) 
becomes considerably simpler in this case. 



III. NUMERICAL IMPLEMENTATION 

In principle, any surface S and any harmonic source 
basis can be used. Point sources, as discussed in Ref. [TJ 
are a simple, although highly inefficient, example. How- 
ever, many common FDTD algorithms (including the one 
we employ in this paper) involve simulation on a dis- 
cretized grid. For these applications, a rectangular sur- 
face S with an expansion basis separately defined on each 
face of S is the simplest. In this case, the field integration 
along each face can be performed to high accuracy and 



converges rapidly. The Fourier cosine series on a discrete 
grid is essentially a discrete cosine transform (DCT), a 
well known discrete orthogonal basis with rapid conver- 
gence properties [23 . This in contrast to discretizing 
some basis such as spherical harmonics that are only ap- 
proximately orthogonal when discretized on a rectangu- 
lar grid. 



A. Two-dimensional systems 

In this section we consider a variant of the piston-like 
configuration of Ref. [8[ shown as the inset to Fig. [3J This 
system consists of two cylindrical rods sandwiched be- 
tween two sidewalls, and is of interest due to the non 
monotonic dependence of the Casimir force between the 
two blocks as the vertical wall separation h/a is varied. 
The case of perfect metallic sidewalls (s(x) = — oo) has 
been solved previously [22 ; here we also treat the case of 
perfect magnetic conductor sidewalls (/jl(x) = — oo) as a 
simple demonstration of method using magnetic materi- 
als. 

While three-dimensional in nature, the system is 
translation-invariant in the ^-direction and involves only 
perfect metallic or magnetic conductors. As discussed 
in Ref. [12] this situation can actually be treated as 
the two-dimensional problem depicted in Fig. [2] using 
a slightly different form for g(—t) in Eq. ([!]) (given 
in the Appendix). The reason we consider the three- 
dimensional case is that we can directly compare the 
results for the case of metallic sidewalls to the high- 
precision scattering calculations of Ref. [22] (which uses 
a specialized exponentially convergent basis for cylin- 
der/plane geometries). 

For this system, the surface S consists of four faces, 
each of which is a line segment of some length L 
parametrized by a single variable x. We employ a co- 
sine basis for our harmonic expansion on each face of S. 
The basis functions for each side are then: 



fn(x) 



fnirx\ 



0,1,... 



(6) 



where L is the length of the edge, and f n (x) = for all 
points x not on that edge of S. These functions, and their 
equivalence to a computation using ^-function sources as 
basis functions, are shown in Fig. [2] 

In the case of our FDTD algorithm, space is discretized 
on a Yee grid [18 , and in most cases x will turn out to 
lie in between two grid points. One can run separate 
simulations in which each edge of S is displaced in the 
appropriate direction so that all of its sources lie on a 
grid point. However, we find that it is sufficient to place 
suitably averaged currents on neighboring grid points, as 
several available FDTD implementations provide features 
to accurately interpolate currents from any location onto 
the grid. 




FIG. 3: Force for the double cylinders of [22] as a func- 
tion of sidewall separation h/a, normalized by the PFA force 
Fpfa = hc((3)d/87ra 3 . Red/blue/black squares show the 
TE/TM/total force in the presence of metallic sidewalls, as 
computed by the FDTD method (squares). The solid lines 
indicate the results from the scattering calculations of [22] , 
showing excellent agreement. Dashed lines indicate the same 
force components, but in the presence of perfect magnetic- 
conductor sidewalls (computed via FDTD). Note that the to- 
tal force is nonmonotonic for electric sidewalls and monotonic 
for magnetic sidewalls. 
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FIG. 4: Tyy. n=2 (t,x) snapshots (blue/ white/red = posi- 
tive/zero/negative) for the n — 2 term in the harmonic cosine 
expansion on the leftmost face of S for the double blocks con- 
figuration of [5] at selected times (in units of a/c). 




The force, as a function of the vertical sidewall sepa- 
ration h/a, and for both TE and TM field components, 
is shown in Fig. [3] and checked against previously known 
results for the case of perfect metallic sidewalls [22 . We 
also show the force (dashed lines) for the case of perfect 
magnetic conductor sidewalls. 

In the case of metallic sidewalls, the force is nonmono- 
tonic in h/a. As explained in Ref. 22, this is due to the 
competition between the TM force, which dominates for 
large h/a but is suppressed for small h/a, and the TE 
force, which has the opposite behavior, explained via the 
method of images for the conducting walls. Switching 
to perfect magnetic conductor sidewalls causes the TM 
force to be enhanced for small h/a and the TE force to be 
suppressed, because the image currents flip sign for mag- 
netic conductors compared to electric conductors. As 
shown in Fig. [3j this results in a monotonic force for this 
case. 

The result of the above calculation is a time-dependent 
field similar to that of Fig. |4j which when manipulated as 
prescribed in the previous section, will yield the Casimir 
force. As in Ref. [TJ our ability to express the force for a 
dissipationless system (perfect-metal blocks in vacuum) 
in terms of the response of an artificial dissipative sys- 
tem (<j ^ 0) means that the fields, such as those shown 
in Fig. |4j rapidly decay away, and hence only a short 
simulation is required for each source term. 

In addition, Fig. [5] shows the convergence of the har- 
monic expansion as a function of n. Asymptotically for 



FIG. 5: Relative contribution of harmonic moment n in the 
cosine basis to the total Casimir force for the double blocks 
configuration (shown in the inset) 

large n, an n -1 / 4 power law is clearly discernible. The 
explanation for this convergence follows readily from the 
geometry of S: the electric field E(x), when viewed as 
a function along S, will have nonzero first derivatives 
at the corners. However, the cosine series used here al- 
ways has a vanishing derivative. This implies that its 
cosine transform components will decay asymptotically 
as n~ 2 [24]. As T E is related to the correlation function 
(E(x)E(x)), their contributions will decay as n -4 . One 
could instead consider a Fourier series denned around the 
whole perimeter of S, but the convergence rate will be 
the same because the derivatives of the fields will be dis- 
continuous around the corners of S. A circular surface 
would have no corners in the continuous case, but on a 
discretized grid would effectively have many corners and 
hence poor convergence with resolution. 



B. Dispersive Materials 

Dispersion in FDTD in general requires fitting an 
actual dispersion to a simple model (eg. a series of 
Lorentzians or Drude peaks). Assuming this has been 
done, these models can then be analytically continued 
onto the complex conductivity contour. 
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FIG. 6: Force per unit length between long silicon waveg- 
uides suspended in air [25], determined by the FDTD method 
(red triangles). Also shown is the analogous one-dimensional 
computation assuming silicon plates of finite thickness (blue 
dashes), and the result F = tt 2 /240<i 4 A for perfect metals 
(assuming plate area A equal to the interaction area of the 
waveguides) . 

As an example of a calculation involving dispersive ma- 
terials, we consider in this section a geometry recently 
used to measure the classical optical force between two 
suspended waveguides [25] , confirming a prediction [26] 
that the sign of the classical force depends on the rel- 
ative phase of modes excited in the two waveguides. 
We now compute the Casimir force in the same geom- 
etry, which consists of two identical silicon waveguides 
in empty space. We model silicon as a dielectric with 
dispersion given by: 

e(u)=e f + g/ "J° 2 (7) 

where coo = 6.6 x 10 15 rad/sec, and eo = 1.035, Sf = 
11.87. This dispersion can be implemented in FDTD 
by the standard technique of auxiliary differential equa- 
tions [18 mapped into the complex-o; plane as explained 
inRef.CD 

The system is translation-invariant in the z direction. 
If it consisted only of perfect conductors, we could use 
the trick of the previous section and compute the force 
in only one 2d simulation. However, dielectrics hybridize 
the two polarizations and require an explicit k z integral, 
as discussed in Ref. [12] Each value of k z corresponds to a 
separate two-dimensional simulation with Bloch-periodic 
boundary conditions. The value of the force for each k z 
is smooth and rapidly decaying, so in general only a few 
k z points are needed. 

To simulate the infinite open space around the waveg- 
uides, it is ideal to have "absorbing boundaries" so that 
waves from sources on S do no reflect back from the 
boundaries. We employ the standard technique of per- 
fectly matched layers (PML), which are a thin layer of ar- 
tificial absorbing material placed adjacent to the bound- 



ary and designed to have nearly zero reflections [T8] . 
The results are shown in red in Fig. [6j We also show 
(in blue) the force obtained using the proximity force 
approximation (PFA) calculations based on the Lifshitz 
formula [27] [28]. For the PFA, we assume two paral- 
lel silicon plates, infinite in both directions perpendic- 
ular to the force and having the same thickness as the 
waveguides in the direction parallel to the force, com- 
puting the PFA contribution from the surface area of the 
waveguide. As expected, at distances smaller than the 
waveguide width, the actual and PFA results are in good 
agreement, while as the waveguide separation increases, 
the PFA becomes more inaccurate. For example, by a 
separation of 300nm, the PFA result is off by 50%. We 
also show for comparison the force for the same surface 
between two perfectly metallic plates, also assuming in- 
finite extent in both transverse directions. 



C. Three dimensions with cylindrical symmetry 

In the case of cylindrical symmetry, we can employ 
a cylindrical surface S and a complex exponential basis 
e im(f) ^ n ^ e (j) direction. For a geometry with cylindrical 
symmetry and a separable source with e %m ^ dependence, 
the resulting fields are also separable with the same (j) de- 
pendence, and the unknowns reduce to a two-dimensional 
(r, z) problem for each m. This results in a substan- 
tial reduction in computational costs compared to a full 
three-dimensional computation. 

Treating the reduced system as a two-dimensional 
space with coordinates (r, z\ the expression for the force 
(as derived in the Appendix) is now: 

Fi = J2 I dt Im k?H)] / <Mx) r ii;n (x,t) (8) 

where the m-dependence has been absorbed into the def- 
inition of T above: 

^ij;n(%it) — ^ij;n,m=0 (x, t) H~2 ^ ^ Re[F^j ;n?m (x, t)] , (9) 

m>0 

and dsj = <isnj(x), ds being a one-dimensional Carte- 
sian line element. As derived in the Appendix, the Ja- 
cobian factor r obtained from converting to cylindrical 
coordinates cancels out, so that the one-dimensional (r- 
independent) measure ds is the appropriate one to use in 
the surface integration. Also, the 2Re[---] comes from 
the fact that the +ra and — m terms are complex conju- 
gates. Although the exponentials e irn ^ are complex, only 
the real part of the field response appears in Eq. ([9|, 
allowing us to use Im[g(—t)} alone in Eq. (J8|. 

Given an e irn ^ dependence in the fields, one can write 
Maxwell's equations in cylindrical coordinates to obtain a 
two-dimensional equation involving only the fields in the 
(r, z) plane. This simplification is incorporated into many 
FDTD solvers, as in the one we currently employ [19] , 
with the computational cell being restricted to the (r, z) 
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FIG. 7: Force as a function of outer sidewall spacing h/a for 
the cylindrically-symmetric piston configuration shown in the 
figure. Both plates are perfect metals, and the forces for both 
perfect metallic and perfect magnetic conductor sidewalls are 
shown. Note that in contrast to Fig. [3] here the force is 
monotonic in h/a for the metallic case and non monotonic for 
the magnetic case. 



plane and m appearing as a parameter. When this is 
the case, the implementation of cylindrical symmetry is 
almost identical to the two-dimensional situation. The 
only difference is that now there is an additional index 
m over which the force must be summed. 

To illustrate the use of this algorithm with cylindrical 
symmetry, we examine the 3d system shown in the inset 
of Fig. [7| This configuration is similar to the configura- 
tion of cylindrical rods of Fig. [3j except that instead of 
translational (z) invariance we instead impose rotational 
((f)) invariance. In this case, the two sidewalls are joined 
to form a cylindrical tube. We examine the force between 
the two blocks as a function of h/a (the h = case has 
been solved analytically [29]). 

Due to the two-dimensional nature of this problem, 
computation time is comparable to that of the two- 
dimensional double block geometry of the previous sec- 
tion. Rough results (at resolution 40, accurate to within 
a few percentage points) can be obtained rapidly on a 
single computer (about 5 minutes running on 8 proces- 
sors) are shown in Fig. [7] for each value of h/a. Only 
indices n, m G {0, 1, 2} are needed for the result to have 
converged to within 1%, after which the error is domi- 
nated by the spatial discretization. PML is used along 
the top and bottom walls of the tube. 

In contrast to the case of two pistons with transla- 
tional symmetry, the force for metallic sidewalls is mono- 
tonic in h/a. Somewhat surprisingly, when the side- 
walls are switched to perfect magnetic conductors the 
force becomes non monotonic again. Although the use 
of perfectly magnetic conductor sidewalls in this exam- 
ple is unphysical, it demonstrates the use of a general- 
purpose algorithm to examine the material-dependence 
of the Casimir force. If we wished to use dispersive 



FIG. 8: The Casimir force between a periodic array of Sil- 
icon waveguides and a Silicon/Silica substrate, as the ar- 
ray/substrate separation is varied. The system is periodic in 
the x-direction and translation-invariant in the ^-direction, 
so the computation involves a set of two-dimensional simula- 
tions. 



and/or anisotropic materials, no additional code would 
be required. 



D. Periodic Boundary Conditions 

Periodic dielectric systems are of interest in many ap- 
plications. The purpose of this section is to demonstrate 
computations involving a periodic array of dispersive sil- 
icon dielectric waveguides above a silica substrate, shown 
in Fig. U 

As discussed in Ref. [12[ the Casimir force for periodic 
systems can be computed as in integral over all Bloch 
wavevectors in the directions of periodicity. Here, as 
there are two directions, x and z, that are periodic (the 
latter being the limit in which the period goes to zero). 
The force is then given by: 



F kz hdk z dk x 



(10) 



where Fk zi k x ls the force computed from one simulation 
of the unit cell using Bloch-periodic boundary conditions 
with wavevector k = (fc^O, k z ). In the present case, the 
unit cell is of period 1/im in the x direction and of zero 
length in the z direction, so the computations are effec- 
tively two-dimensional (although they must be integrated 
over k z ). 

We use the dispersive model of Eq. ([7]) for silicon, while 
for silica we use [8] 



e(u) = 1 • 



3 

E 

3 



Cjoj 



(11) 



where (C 1 ,C 2 ,C 3 ) = (0.829,0.095,1-098) and 
(lji,lj 2 ,W3) = (0.867,1.508,203.4) x 10 14 (rad/sec). 
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FIG. 9: Three-dimensional configuration showing stable levi- 
tation. At the equilibrium point, the force of gravity counters 
the Casimir force, while the Casimir force from the walls of 
the spherical indentation confine the sphere laterally 
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FIG. 10: Total (Casimir + gravity) vertical (z) force on the 
silica sphere (depicted in the inset) as the height h of the 
sphere's surface above the indentation surface is varied. The 
point of vertical equilibrium occurs at h ~ 450 nm. 



E. Full 3d computations 

As a final demonstration, we compute the Casimir 
force for a fully three-dimensional system, without the 
use of special symmetries. The system used is depicted 
in Fig. [§} 

This setup demonstrates stable levitation with the aid 
of a fluid medium, which has been explored previously 
in Ref. 9. With this example, we present a setup simi- 
lar to that used previously to measure repulsive Casimir 
forces [6], with the hope that this system may be exper- 
imentally feasible. 

A silica sphere sits atop a perfect metal plane which has 
a spherical indentation in it. The sphere is immersed in 
bromobenzene. As the system satisfies £ sp here < £fluid < 
£piane, the sphere feels a repulsive Casimir force up- 
wards [6 . This is balanced by the downward force of 
gravity, which confines the sphere vertically. In addition, 
the Casimir repulsion from the sides of the spherical in- 
dentation confine the sphere in the lateral direction. The 
radius of the sphere is 1/im, and the circular indenta- 
tion in the metal is formed from a circle of radius 2/im, 
with a center 1/im above the plane. For computational 
simplicity, in this model we neglect dispersion and use 
the zero-frequency values for the dielectrics, as the basic 
effect does not depend upon the dispersion (the precise 
values for the equilibrium separations will be changed 
with dispersive materials). These are e = 2.02 for silica 
and e = 4.30. 

An efficient strategy to determine the stable point is 
to first calculate the force on the glass sphere when its 
axis is aligned with the symmetry axis of the indenta- 
tion. This configuration is cylindrically-symmetric and 
can be efficiently computed as in the previous section. 
Results for a specific configuration, with a sphere radius 
of 500 nm and an indentation radius of 1 /im, are shown 
in Fig. fT0] 
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FIG. 11: Casimir restoring force on the sphere as a function of 
lateral displacement dx, when the vertical position is fixed at 
h = 450 nm, the height at which gravity balances the Casimir 
force 



The force of gravity is balanced against the Casimir 
force at a height of h = 450 nm. To determine the 
strength of lateral confinement, we perform a fully three- 
dimensional computation in which the center of the 
sphere is displaced laterally from equilibrium by a dis- 
tance dx (the vertical position is held fixed at the equilib- 
rium value h = 450nm). The results are shown in Fig. 11 
It is seen that over a fairly wide range (|Ax| < 100 nm) 



the linear term is a good approximation to the force, 
whereas for larger displacements the Casimir force be- 
gins to increase more rapidly. Of course, at these larger 
separations the vertical force is no longer zero, due to the 
curvature of the indentation, and so must be re-computed 
as well. 

The fully three-dimensional computations are rather 
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large, and require roughly a hundred CPU hours per 
force point. However, these Casimir calculations paral- 
lelize very easily - every source term, polarization, and k- 
point can be computed in parallel, and individual FDTD 
calculations can be parallelized in our existing software - 
so we can compute each force point in under an hour on 
a supercomputer (with 1000+ processors). In contrast, 
the 2d and cylindrical calculations require tens of minutes 
per force point. We believe that this method is usable 
in situations involving complex three-dimensional mate- 
rials (e.g, periodic systems or systems with anisotropic 
materials) . 



IV. CONCLUDING REMARKS 

We have demonstrated a practical implementation of 
a general FDTD method for computing Casimir forces 
via a harmonic expansion in source currents. The utility 
of such a method is that many different systems (disper- 
sive, anisotropic, periodic boundary conditions) can all 
be simulated with the same algorithm. 

In practice, the harmonic expansion converges rapidly 
with higher harmonic moments, making the overall com- 
putation complexity of the FDTD method 0(iV 1+1 / d ) 
for N grid points and d spatial dimensions. This arises 
from the O(N) number of computations needed for one 
FDTD time step, while the time increment used will 
vary inversely with the spatial resolution [18 , leading 
to 0(N 1 ^ d ) time steps per simulation. In addition, there 
is a constant factor proportional to the number of terms 
retained in the harmonic expansion, as an independent 
simulation is required for each term. For comparison, 
without a harmonic expansion one would have to run a 
separate simulation for each point on S. In that case, 
there would be 0{N ( ^ d ~ 1 ^ d ) points, leading to an overall 
computational cost of 0(N 2 ) pQ. 

We do not claim that this is the most efficient tech- 
nique for computing Casimir forces, as there are other 
works that have also demonstrated very efficient meth- 
ods capable of handling arbitrary three-dimensional ge- 
ometries, such as a recently-developed boundary-element 
method [13]. However, these integral-equation methods 
and their implementations must be substantially revised 
when new types of materials or boundary conditions are 
desired that change the underlying Green's function (e.g., 
going from metals to dielectrics, periodic boundary con- 
ditions, or isotropic to anisotropic materials), whereas 
very general FDTD codes, requiring no modifications, 
are available off-the-shelf. 
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VI. APPENDIX 

A. Simplified computation of g(t) 

In Ref. 1 we introduced a geometry-independent func- 
tion g(t), which resulted from the Fourier transform of a 
certain function of frequency, termed which is given 
by ® : 



9(0 



(12) 



Once g(t) is known, it can be integrated against the fields 
in time, allowing one to compute a decaying time-series 
which will, when integrated over time, yield the correct 
Casimir force. 

g(£) has the behavior that it diverges in the high- 
frequency limit. For large £, g(£) has the form: 



0(0 -01(0^6(0 ■ 



•0-6(0 as^oo (13) 



Viewing #i(0 as a function, we could only compute its 
Fourier transform g±(t) by introducing a cutoff in the 
frequency integral at the Nyquist frequency, since the 
time signal is only defined up to a finite sampling rate 
and the integral of a divergent function may appear to 
be undefined in the limit of no cutoff. 

Applying this procedure to compute g(—t) yields a 
time series that has strong oscillations at the Nyquist 
frequency. The amplitude of these oscillations can be 
quite high, increasing the time needed to obtain conver- 
gence and also making any physical interpretation of the 
time series more difficult. 

These oscillations are entirely due to the high- 
frequency behavior of where ~ <7i(0- However, 
g(t) and g(£) only appear when when they are being in- 
tegrated against smooth, rapidly decaying field functions 
r(x, t) or r(x, £). In this case, g can be viewed as a 
tempered distribution (such as the ^-function) [30]. Al- 
though g(£) diverges for large £, this divergence is only a 
power law, so it is a tempered distribution and its Fourier 
transform is well-defined without any truncation. In par- 
ticular, the Fourier transform of #i(£) is given by: 



2tt \t 2 



(14) 



Adding and subtracting the term gi(£) from the 
remaining term decays to zero for large £ and can be 
Fourier transformed numerically without the use of a 
high-frequency cutoff, allowing g(—t) to be computed as 
the sum of gi(t) plus the Fourier transform of a well- 
behaved function. This results in a much smoother g(—t) 
which will give the same final force as the g(—t) used 
in Ref. 1, but will also have a much more well-behaved 
time dependence. 

In Fig. [12] we plot the convergence of the force as a 
function of time for the same system using the g(—t) 
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FIG. 12: Plot of the force error (force after a finite time in- 
tegration vs. the force after a very long run time) for g(t) 
determined from a numerical transform as in Ref . [T] and from 
the analytic transform of the high-frequency components. In- 
set: lm[g{— £)] obtained without a cutoff, in which the high- 
frequency divergence is integrated analytically. Compare with 
Fig. 1 of Ref. [1] 



obtained by use of a high-frequency cutoff and for one 
in which #i(£) is transformed analytically and the re- 
mainder is transformed without a cutoff. The inset plots 
Im<?(— t) obtained without using a cutoff (since the real 
part is not used in this paper) for a = 10. If a complex 
harmonic basis is used, one must take care to use the full 
g(t) and not only its imaginary part. 



1. Further simplification 

In addition to the treatment of the high-frequency di- 
vergence in the previous section, we find it convenient 
to also Fourier transform the low-frequency singularity 
of g(£) analytically. As discussed in Ref. TJ the low- 
frequency limit of g(£) is given by: 



g(0 - 92(0 



Via 3 / 2 
2 £V2 



9(£) as £ -> (15) 



The Fourier transform of §2 (£)> viewed as a distribu- 
tion, is: 



92(~t) 



T 3/2 



4v^F tV2 



(16) 



After removing both the high- and low-frequency di- 
vergences of we perform a numerical Fourier trans- 
form on the function 5g(£) = g (£) — g\ (£) — #2 (0 > which is 
well-behaved in both the high- and low- frequency limits. 

In the present text we are only concerned with real 
sources, in which case all fields r(x, t) are real and only 
the imaginary part of g(—t) contributes to the force 
in Equation ([TJ. The imaginary part of g(—t) is then: 



lm[g(-t)} = lm(Sg(-t)) 



1 

2^ 



1 a 



3/2 



4v^ t 1 / 2 
(17) 



2. Perfect conductors and z-invariance 

As discussed in Ref. [12j the stress-tensor frequency 
integral for a three-dimensional ^-invariant system in- 
volving only vacuum and perfect metallic conductors is 
identical in value to the integral of the stress tensor for 
the associated two-dimensional system (corresponding to 
taking a z = crossection), with an extra factor of ico/2 
in the frequency integrand. In the time domain, this cor- 
responds to solving the two-dimensional system with a 
new g(—t). 

In this case the Fourier transform can be performed 
analytically. The result is: 



lm[g(-t)} 



-(- 



3ct 
2fi 



a 
2t 



(18) 



B. Harmonic expansion in cylindrical coordinates 

The extension of the above derivation to three dimen- 
sions and non-Cartesian coordinate systems is straight- 
forward, as the only difference is in the representation of 
the (5-function. Because the case of rotational invariance 
presents some simplification, we will explicitly present 
the result for this case below. 

For cylindrical symmetry, we work in cylindrical coor- 
dinates (r, (/>, z) and choose a surface S that is also rota- 
tionally invariant about the z-axis. S is then a surface of 
revolution, consisting of the rotation of a parametrized 
curve (r(s), <fi = 0, z(s)) about the z axis. The most prac- 
tical harmonic expansion basis consists of functions of the 
form fn(x)e irn ^. Given a (j) dependence, many FDTD 
solvers will solve a modified set of Maxwell's Equations 
involving only the (r, z) coordinates. In this case, for each 
m the problem is reduced to a two-dimensional problem 
where both sources and fields are specified only in the 
(r, z)-plane. 

Once the fields are determined in the (r, z)-plane, the 
force contribution for each m is given by: 



2tt 



[ 2l> V / C Z S (x')r(x')e im *'(55(x-x')rg ;m (r;x,x') 
Jo Js 



) / dsj(x) r(x)e 
Js 



imcj) 



(19) 



where the values of x range over the full three- 
dimensional (r, </>, z) system. Here we introduce the 
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Cartesian line element ds along the one-dimensional sur- 
face S in anticipation of the cancellation of the Jacobian 
factor r(x) from the integration over S. We have explic- 
itly written only the contribution for T E , the contribution 
for being identical in form. 

In cylindrical coordinates, the representation of the 5- 
function is: 



5(x - x') = 



1 



27rr(x) 



5(0 - <l>')8{r - r')5(z - z') (20) 



For simplicity, assume that S consists entirely of z = 
constant and r = constant surfaces (the more general 
case follows by an analogous derivation). In these cases, 
the surface 5-function 5s is given by: 



*s(x-x') 
Ss(x-x') 



27rr(x) 



27rr(x) 



- 4>)5{r — r'), z = constant 
S((j) — 4>')5{z — z'\ r = constant 



In either case, we see that upon substitution of either 
form of 5s into Eq. (19), we obtain a cancellation with 



the first r(x) factor. Now, one picks an appropriate de- 
composition of 5s into functions f n (a choice of r = const 
or z = const merely implies that the f n will either be 
functions of z, or r, respectively). We denote either case 
as / n (x), with the r and z dependence implicit. 

We now consider the contribution for each value of n. 
The integral over x' is: 



the case of cylindrical symmetry, this field must have a 
(j) dependence of the form e im ^: 



(21) 



This factor of e im ^ cancels with the remaining e _zm ^. 
The integral over (j) then produces a factor of 2ir that 
cancels the one introduced by 5s- After removing these 
factors, the problem is reduced to one of integrating the 
field responses entirely in the (r, z) plane. The contribu- 
tion for each n and m is then: 



L 



dsj(x) /„(x)rg.„ m (t,r, z) 



(22) 



If one chooses the / n (x) to be real- valued, the contri- 
butions for +ra and — m are related by complex conjuga- 
tion. The sum over m can then be rewritten as the real 
part of a sum over only non negative values of m. The 
final result for the force from the electric field terms is 
then: 

Fi = f dthn\g(-t)]J2 [ ^(r,z)/ n (r,z)rg ;n (t,r,z) 
Jo n Js 

(23) 

where the m-dependence has been absorbed into the def- 
inition of Tij. n as follows: 



r 2ir 



^0 Js 



d S (x')r(x')r^ ;nm (t,x,x')/„(x')e 



As noted in the text, r^. nm (t,x) is simply the field 
measured in the FDTD simulation due to a three- 
dimensional current source of the form / n (x)e zm( ^. In 



^ij;n(^5 T i Z ) — ^ij;n,m=o(^i z) -\- 2 ^ ] Re[r^ ;nm (t, T, z)\ 

m>0 

(24) 

We have also explicitly included the dependence on r 
and z to emphasize that the integrals are confined to 
the two-dimensional (r, z) plane. The force receives an 
analogous contribution from the magnetic-field terms. 
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